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Abstract 



We analyze in detail the error that arises from the linearization in linearized augmented-plane-wave (LAPW) basis functions around 
predetermined energies Ei and show that it can lead to undesirable dependences of the calculated results on method-inherent 
parameters such as energy parameters Ei and mufhn-tin sphere radii. To overcome these dependences, we evaluate approaches 
,that eliminate the linearization error systematically by adding local orbitals (LOs) to the basis set. We consider two kinds of 
■LOs: (i) constructed from solutions M/(r, E) to the scalar-relativistic approximation of the radial Dirac equation with E > Ei and (ii) 
constructed from second energy derivatives d^ui{r, E)/dE^ at £ = Ej. We find that the latter eliminates the error most efficiently and 
yields the density functional answer to many electronic and materials properties with very high precision. Finally, we demonstrate 
that the convergence behavior of the so constructed LAPW+LO basis is superior to that of the related APWh-Io approach in cases 
where the linearization error is large. 
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1. Introduction 

In the past decades, material simulations have become an in- 
valuable approach in condensed matter physics and materials 
|science. Ever increasing computer power as well as theoretical 
and methodological progress in the description of materials are 
the main incentives for more and more accurate calculations on 
more and more complex materials. Within the wide range of 
'theoretical approaches, density functional theory (DFT) yj] is 
the method of choice for the calculation of electronic ground- 
state properties of materials. 

' Practical realizations of DFT almost invariably rely on the 
Kohn-Sham (KS) formalism 12], which employs an auxiliary 
system of noninteracting electrons whose number density co- 
incides with that of the real interacting system. Most codes 
make use of a set of basis functions to represent the quantum 
mechanical wave functions of these noninteracting electrons, 
■which enables a formulation of the underlying differential KS 
equation as a generalized eigenvalue problem. In recent years 
■we witnessed a trend toward the investigation of solids of in- 
creasing electronic, chemical, and structural complexity, solids 
that exhibit narrow electronic bands, large band gaps, electrons 
that contribute to physisorption and chemisorption, which are 
in turn described by more sophisticated exchange and corre- 
lation functionals, e.g., hybrid functionals, the exact-exchange 
functional in the optimized-effective-potential method, van der 
Waals functionals, or correlation functionals based on the ran- 
dom phase approximation (RPA), just to name a few. The 
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description of the electronic structure and the application of 
the new functionals provide new challenges to the efficiency 
and ability of the basis set to precisely represent the density- 
functional answer The aim of this paper is to evaluate and im- 
prove the LAPW basis for this purpose. 

The simplest basis set for systems with periodic boundary 
conditions is certainly the plane-wave basis. The accuracy of 
which is controlled by a single convergence parameter, the mo- 
mentum cutoff radius G,nax- However, the rapid variations of 
the wave functions close to the atomic nuclei cannot be resolved 
in practice with this basis, and one has to resort to pseudopo- 
tentials and pseudized wave functions within a certain distance 
from the atomic nuclei Istl- The core electrons are then incorpo- 
rated into the pseudopotentials, and only the valence electrons 
are treated explicitly. 

The pseudopotential approximation effectively restricts the 
range of materials that can be examined. Compounds contain- 
ing 4/ and 5 / elements, and transition-metals as well as ox- 
ides, nitrides, and carbides cannot be treated efficiently within 
this approach. Among the all-electron approaches that de- 
scribe core and valence electrons on an equal footing (Gaus- 
sian functions HI], the projector augmented-wave Jsl and lin- 
earized muffin-tin orbital method Q |3> Si to name a few), 
the full-potential linearized augmented-plane-wave (FLAPW) 
method |0,|^[l3l provides one of the most precise basis sets for 
all-electron calculations. It allows for studying the electronic 
structure of a large variety of materials, including open systems 
with low symmetry and compounds of any chemical composi- 
tion. 

The FLAPW method is based on a partitioning of space into 
non-overlapping spheres centered at the atomic nuclei, the so- 
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called muffin-tin (MT) spheres, and the interstitial region (IR). 
The core states are completely confined within the spheres, 
which allows to treat them as localized states in a spherically 
symmetric atomic potential. For the valence electrons, on the 
other hand, the APW basis functions are defined piecewise: 
plane waves of all reciprocal lattice vectors up to the maxi- 
mal momentum Gmax in the IR, which are augmented by radial 
functions in the MT spheres that are solutions of the scalar- 
relativistic approximation to the Dirac equation for the spher- 
ically averaged effective potential. Employing the linearized 
APW (LAPW) basis set, the energy-dependent radial functions 
are approximated by energy-independent functions evaluated at 
predetermined energy parameters £/. The functions are linearly 
combined so as to match to the plane waves in value and slope 
at the MT sphere boundaries. In the conventional LAPW ba- 
sis there are two radial functions per angular momentum quan- 
tum number /, the solution M;(r, £/) and its energy derivative 
ui{r,Ei) - dui{r, E)/dE\E=Er The inclusion of the derivative 
provides variational freedom around the energy parameters . 

This conventional LAPW basis set is a very accurate one for a 
wide range of materials. In comparison to the pure plane-wave 
basis, the LAPW basis requires far less basis functions, while 
still being able to provide an all-electron description. However, 
no matter how large the momentum cutoff radius is chosen, the 
flexibility of the basis in the MT spheres is restricted to the two 
radial functions M;(r, £/) and M/(r, £;) per I quantum number and 
the energy-dependent radial function is approximated by these 
energy-independent functions. This gives rise to a lineariza- 
tion error that can notably affect the accuracy, for example, for 
materials with large bandwidths, large bandgaps, or if states 
are considered that are energetically far away from the energy 
parameters OI ll Il2ll . It is obvious that the linearization error 



depends on two sets of method-inherent parameters, (i) the en- 
ergy parameters Ei and (ii) the MT radii, Rmt, as they define 
the region of space in which the wave function is represented 
by El) and M;(r, £/). None of the parameters can be cho- 
sen on the basis of a variational principle. Optimally, the final 
results should not depend on Ei and Rmt as long as they are 
chosen within a reasonable range of values. 

Over time, several approaches have been proposed to reduce 
the linearization error. The first approach that we mention here 
is the quadratic APW (QAPW) method OQ- In this method 
the MT augmentation is extended by including the second-order 
energy derivative M/(r, £/) and employing an algebraic relation 
for the matching coefficients of M/(r, £/) and M;(r, £/). 

Next, Singh 1 15 . T^l investigated how to deal with semicore 
states. These are high-lying core states that are not completely 
confined within the MT spheres and therefore cannot be treated 
separately from the valence states. He introduced the radial 
functions ui{r, Ef''') that solve the scalar-relativistic Dirac equa- 
tion with energy parameters Ef^ <sc Ei that correspond to the 
energy of the respective semicore state. He introduced and 
compared the inclusion of these functions by either enforcing 
continuity of the basis-function curvature as additional match- 
ing condition at the MT boundaries or constructing additional 
local orbitals (LOs) with M/(r, £/), M/(r, £;), and ui{r,E^'~) that 
are nonzero only in the MT spheres. He found that the first ap- 



proach is effective in describing the semicore states, on the ex- 
pense of creating stiffer basis functions across the MT boundary 
imposed by the additional matching conditions. As a result, the 
basis set becomes less flexible and larger basis sets are required 
to achieve the same accuracy. On the other hand, the extension 
of the basis with LOs avoids this problem. 

The extended LAPW (ELAPW) basis developed by 



Krasovskii et al. IIITl II8L II9I1 is another approach amending 



the conventional LAPW basis set by LOs to improve the de- 
scription of the electronic structure over a wide energy win- 
dow. Here, pairs of LOs are constructed from M/(r, E^'^) and 
M/(r, E^^), respectively, at energy parameters £'J"° that are typi- 
cally placed at energies above the Fermi energy. However, the 
authors did not give a prescription of how to choose them in an 
optimal way. Friedrich et al. |fl2ll considered LOs that are con- 
structed from second-order energy derivatives M/(r, £;) taken at 
the valence energy parameters £/ to improve the description of 
high-lying states for calculations. Betzinger ef a/. 12011 em- 
ployed LOs defined at higher energies to refine the susceptibil- 
ity matrix in optimized-effective-potential calculations. 

Sjostedt et al. 112 111 used LOs to achieve a linearization of the 
basis set, known as the APW+lo method, in which the con- 
dition of continuous slope across the MT sphere boundary is 
dropped for the basis functions. The functions H/(r, Ei) are not 
used anymore as augmentation functions but instead are in- 
cluded in LOs. This makes the basis set more flexible when 
compared with the FLAPW method. Thus, the same accuracy 
can be achieved with less basis functions and this saves com- 
putation time. We will show later that the inclusion of LOs, 
in particular second-derivative LOs, into the LAPW basis can 
yield a similar gain in flexibility. With respect to accuracy the 
resulting basis is even superior to APW+lo in cases where the 
linearization error is large. 

In this article we evaluate the conventional LAPW basis set 
for a set of systems exhibiting large muffin-tin radii and large 
band gaps. We have chosen fee Ce, KCl in the rock-salt struc- 
ture, fee Ar and bcc V as test systems. We examine the lin- 
earization error and find that physical quantities such as the to- 
tal energies, lattice constants, or band structures are susceptible 
to slight variations of method-inherent parameters such as MT 
radii and energy parameters. We apply two types of local or- 
bitals to extend the LAPW basis set, (i) one based on functions 
of higher energy derivative (LAPW+HDLO) and (ii) one based 
on energy parameters at higher energy (LAPW-i-HELO). We 
show that these LOs are very effective in reducing or, within the 
relevant degree of accuracy, even in eliminating the lineariza- 
tion error. Here, the second-derivative HDLOs are particularly 
efficient. 

The article continues with a brief introduction to the LAPW 
basis and the definition of the LO extensions in section |2l In 
section [3j we introduce the test set of materials for which our 
evaluations are carried out. We continue in section |4] with a 
discussion of the Unearization error for the conventional and 
the LO extended LAPW basis sets. In detail, in section 143] we 
evaluate the quality of the various basis sets in their ability to 
represent Kohn-Sham states inside the MT sphere, discuss de- 
pendences of computed physical properties on method-inherent 
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parameters in section 14.21 and investigate the effect of the lin- 
earization error on unoccupied states, in particular, the KS band 
gap in section 143] Finally, we discuss in section[5]the effect of 
the LOs on the convergence behavior of the total energy with 
respect to the basis set size for the different basis sets and also 
in comparison with the APW+lo approach, before we conclude 
in section|6l 

2. LAPW basis and LO extensions 

In the KS formalism of DFT Q] one considers a fictitious 
system of noninteracting electrons moving under the influence 
of an effective potential V"'*'^(r). For lattice periodic potentials 
the electronic quantum mechanical wave functions i^kn(r) with 
the Bloch vector k and the band index n are thus solutions of 



--A + y'='^(r) 
2 



(1) 



where the spin index has been suppressed. The effective poten- 
tial is defined such that the KS system exhibits the same number 
density as that of the real interacting system. It comprises the 
electrostatic potential created by the charged particles - elec- 
trons and nuclei - as well as an exchange-correlation potential, 
for which we use here the Perdew-Zunger parametrization of 
the local-density approximation ll22ll . For simplicity, we present 
in Eq. ([TJ and the following the non-relativistic equations, while 
in practice the scalar- or fully relativistic Dirac equation is em- 
ployed. 

Close to the atomic nuclei the effective potential is predomi- 
nantly spherical, which allows the core states to be determined 
efficiently from the fully relativistic radial Dirac equation. The 
core states fall off quickly toward the MT sphere boundary. 



where they are practically zero. It can be shown llI2[ Il6l that 



this property guarantees that the core states are orthogonal to 
the radial functions used for the construction of the LAPW basis 
for the valence electrons, which we introduce in the following. 
The valence states are represented by the set of basis functions 



0kG(r) 



_L „/(k+G)r 



for r e IR 



2 R'^^ir,, EiMh) for r € MT„ 



(2) 



where the G are reciprocal lattice vectors, Q. is the unit-cell vol- 
ume, L = (/, m) comprises the angular momentum and the mag- 
netic quantum number, r„ is the position vector relative to the 
MT sphere center R'' of atom a, and Yi{r) are the spherical har- 
monics. In practice, one employs the cutoff parameter Gmax that 
controls the number of basis functions through |k + G| < Gmax 
and implies the cutoff parameter /maxo- with / < /maxo- through 
the rule of thumb Z^axa = GmaxRuTa Ell- 

For each angular momentum quantum number I, the radial 
functions 

R^i^ir,, El,) = fl^^Mto(r„ + b^^ii,,(r,, (3) 

are linear combinations of the normalized solution to 

f 1 d- l{l + 1) 



+ Vf(r) \ r uUr, E) ^ E r u,„{r, E) (4) 



2 dr^- 2r2 



for r e MTq, and its energy derivative uia{r, Eia), where E - Eia 
and V^ir) is the spherical part of V^^ir) inside of MT„ and ex- 
pressed in Hartree units (1 Htr = 2 Ry a: 27.2 eV). The matching 
coefficients and are determined such that the 0kG(r) are 
continuous in value and slope at the MT sphere boundaries. 

The energy is a parameter that is fixed for each iteration 
of a self-consistent-field run. If it happens to be identical to the 
KS eigenvalue et,,, the function uia(r,Eia)YL{r) already solves 
Eq. ([TJ pointwise in the MT sphere by construction (provided 
that we restrict ourselves to the spherical part of the effective 
potential, which is indeed much larger than the nonspherical 
terms). In this sense, the inclusion of the energy derivative 
uiair, Eia) makes it possible that states with energies different 
from but close to Eia can also be described accurately, in accor- 
dance with the Taylor expansion to linear order 



Ula(_r, E) a uiair, Eia) + {E - Eia)uia{r, Eia) ■ 



(5) 



The energy parameters Eia should be chosen as close as pos- 
sible to the corresponding band energies. While the bands of a 
crystal exhibit a strong dependence on the Bloch vector k and 
the band index n, the energy parameters only depend on the 
atom and the angular momentum /. Thus, differences between 
the band energies and the energy parameters are unavoidable. 

In practice, several methods are in use to choose the en- 
ergy parameters automatically in each iteration of the self- 
consistent-field cycle. One example for such a method solves 
for each sphere an atomic Hamiltonian employing the spherical 
part of the effective KS potential in the MT sphere with a con- 
fining potential outside. The energy parameters of the valence 
states are then set to the corresponding atomic eigenenergies. 
The so-chosen energy parameter, that we denote as atomic en- 
ergy parameter, is motivated by the fact that the bands in a solid 
form out of the atomic states. However, this approach may yield 
energy parameters above the Fermi energy since it does not take 
the occupation of states into account. Energy parameters below 
the Fermi energy are obtained by another method that sets them 
at the energy center of mass of the /-resolved partial density 
of the occupied valence states. The disadvantage of the latter 
method is that systems with narrow bands at the Fermi energy, 
where changes in the occupation of bands can occur easily be- 
tween two subsequent self-consistent-field iterations, exhibit a 
stronger variation of the energy parameters, which at the end 
may reduce the speed to self-consistency or even make the self- 
consistency process less stable than the former choice of £/, 
which is less sensitive to these band reorderings. Both methods 
also depend on the chosen MT radiiQ 

Of course, the more the KS eigenvalues differ from the en- 
ergy parameter, the less adequate the basis becomes for the cor- 
responding KS eigenfunctions, which we refer to as the lin- 
earization error From Eq. (|5) a solution seems obvious: one 
simply adds the second-order derivative uia{r, Eia) to Eq. Q. 
However, this would require an additional expansion coefficient 



For completeness we note that we employ the atomic energy parameters 
for El if not stated otherwise. The qualitative results do not depend on the 
particular choice of the parameters. 
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and, as a consequence, an additional matching condition at the 
MT boundaries. It was shown OlSll that this leads to a less flex- 
ible (slower convergent) basis set. The LO construction u5\\ 
avoids this problem. 

LOs are additional basis functions 



(6) 



that are completely confined to a MT sphere. There are 2/ + 1 
LOs per / quantum number Their radial part is a linear combi- 
nation 



where the coefficients are determined by enforcing normaliza- 
tion J r^R^^(ra)^dra = 1 as well as zero value and slope at the 
MT boundary, implying that the LOs are zero in the IR. The 
third radial function M^°(r„) can be chosen rather arbitrarily. 
For a semicore state it should be the solution of Eq. ([4} with 
the eigenenergy of the semicore state as the energy parameter 
for the given angular momentum / jlSll . The conduction-band 
spectrum may be improved by adding LOs with ener gy p aram- 
eters chosen in the corresponding energy range 

MM- The 

third radial function can also be taken to be the second-order en- 
ergy derivative iiiair, E ig), which would correspond to the next 
higher term in Eq. (|5]l lll2ll . 

We consider here two kinds of LOs: (i) Higher-derivative 
LOs (HDLOs), where the third radial function u^^{r) is set 
equal to iiiair, Eia) (here we only consider the second deriva- 
tive) and (ii) Higher-energy LOs (HELOs), where Mf°(r) = 



uia{r,E\°) with E]'^ > Ei^. For the HELOs we determine E\° 
by the condition that the logarithmic derivative 



,'LO, 



D,a(E) 



(r,E) 



u\-°(r,E) 



"la 



(8) 



is equal to the constant -(/+ 1), corresponding to the evanescent 
solution of Eq. (|4]i with V^** = £ = 0. Because of the form of 
Dia{E) there are infinitely many solutions that are all orthogonal 
to each other and each being related to a particular principal 
quantum number This allows a systematic extension of the 
basis. Figure [Ha) displays the logarithmic derivatives for the 
d channel of fee cerium as an example, where the intersection 
with Di = -3 gives the set of energy parameters E^^ for Ce. 
For reference, we plot in Fig.[Tlb) the radial function m/=2('", E) 
for different energies E covering a range of L8 Htr. The pole 
of Di at around E - 0.4 Htr in Fig. [Tta) marks the transition 
from 5d to 6d orbitals and is related to the zero value of the 
corresponding radial function at the sphere boundary, the node 
at around E = 1.3 Htr in Fig. [Tta) corresponds to a vanishing 
slope in (b). 

Fig.[TJb) displays the energy dependence of the Ce 5d wave 
function, M/(r, E), inside the MT sphere. In the vicinity of the 
atomic nucleus up to a radius of about 1 .5 ao, ui{r, E) is not very 
sensitive to the particular choice of energy E, but it is deter- 
mined by the singularity of the effective potential and the angu- 
lar momentum barrier This changes towards the MT boundary, 
where chemical bonding determines the details of the potential. 
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Figure 1: (a) Logarithmic derivative Di{E) for the d channel of fee Ce as a 
function of the energy relative to the 5;^ energy parameter Esj. The energy 
parameters for the two additional sets of LOs in the LAPW+HELOxl and 
LAPW+HELOx2 basis sets are marked, (b) Solutions of Eq. @ obtained for 
different energies E. The functions are not normalized but scaled to coincide at 
the first maximum. The energy dependence increases towards the MT sphere 
boundary at 3.14 ciq- 



and a strong energy dependence of the basis function is dis- 
played. As a consequence, one can expect that the choice of the 
MT radius may have a significant influence on the linearization 
error With increasing MT radius we expect an increase of the 
error On the other hand a larger MT radius is advantageous 
since it implies a smaller number of required basis functions or 
GmsLx, respectively. We mention in passing that a reduction of 
the MT radius broadens the branches of the logarithmic deriva- 
tive. This also shifts the HELO energies upwards. 

In this work, we investigate and compare HDLO and HELO 
extensions to the LAPW basis. We work with two different 
sets of HELOs that include one and two sets of LOs of suc- 
cessive principle numbers, denoted by HELOx 1 and HELOx2, 
respectively. The LAPW+HDLOxl and LAPW-i-HELOxl ba- 
sis sets comprise 16 additional LOs per atom with / = 0, 1, 2, 3, 
while the LAPW-i-HELOx2 basis contains 32 additional LOs 
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DOS (states / eV) parameters 

Figure 2: Density of states and band structure of fee Ce. Tlie left panel shows 
the band structure along the high-symmetry lines K-F-L. The center panel 
shows the associated total density of states, as well as its projection onto the 
MT sphere and onto the respective .v, p, d, and / channels. The right panel 
displays the atomic energy parameters. 

per atom. 

3. Investigated materials and calculational parameters 

The quality of the different basis sets is analyzed by studying 
the linearization error and its consequences in terms of physical 
quantities for a test set of representative materials, i.e. materi- 
als for which a significant linearization error can be expected. 
This includes materials with large MT radii (around 3 ao), large 
band gaps (between 5-10 eV), or cases where the choice of en- 
ergy parameters deviates significantly from the electronic state 
to be described. In total we have selected fee Ce, KCl in the 
rock-salt structure, fee Ar, and bcc V. Typical cut-off parame- 
ters such as the number of k points in the irreducible wedge of 
the Brillouin zone (IBZ) or the number of basis functions per 
atom controlled by the product Gmax^MX are chosen such that 
the quantities discussed in the paper are converged with respect 
to these parameters. 

The first material to be taken under scrutiny is fee cerium. 
Its choice is largely motivated by previous convergence stud- 
ies of APW-type basis sets documented in literature [i2l]| . The 
MT radius is typically chosen to be larger than 3 qq (uq is the 
Bohr radius). In the standard configuration we perform non- 
magnetic calculations with the experimental lattice constant of 
flexp - 9.05 flo, a MT radius of Rmt - 3.14 ao, GmaxRwn - 13 
corresponding to 222 LAPW basis functions per atom, and 
Imax = 12. As we will see in section |5] this will lead to abso- 
lute convergence of the total energy of about 14.4 meV for the 
conventional LAPW basis, of about 3.3 meV for the HELOxl 
extended basis, and less than 1 meV for the other LO extended 
basis sets. Energy differences are converged to an accuracy that 
is one order of magnitude higher than required in the test cal- 
culations. We describe the 4/, 5d, 6s, and 6p valence electrons 
with the LAPW basis and also include the 5 s and 5/:' semicore 
states by additional LOs in the valence basis. Moreover, we use 
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Figure 3: Density of States and atomic energy parameters for KCl. 



182 k points in the IBZ and atomic energy parameters. Fig- 
ure |2] shows the density of states (DOS) and the band structure 
along the high-symmetry lines K-F-L together with the energy 
parameters employing the conventional LAPW basis. The ob- 
tained energy parameters all lie within their associated bands 
and within a distance of only a few eV to the occupied valence 
states that cover an energy interval of 3 eV below the Fermi 
level. 

The second material, KCl in the rock-salt structure, features 
a large KS band gap of about 5 eV. To describe this gap ac- 
curately, electronic states over a broad energy range have to be 
covered. In the standard configuration we use the same radii of 
Rmt - 2.8 ao for both types of atoms and perform the calcula- 
tions with the experimental lattice constant of Oexp - 11.89 oq. 
Furthermore, we use Gmax^MX = 13 (~355 LAPW basis func- 
tions per atom), Zniax =12, and 60 k points in the IBZ. The 3 s 
and 3p states in K are described by means of semicore LOs in 
the valence-band window. Figure [3] shows the density of states 
together with the energy parameters. The band structure is in- 
vestigated in detail in section l43] and is shown in figure[T3] Be- 
yond the energy range shown in the figures, the energy parame- 
ter for the 4/ states in K is positioned at E4 / k -Ep - 15.63 eV, 
the parameter for the CI 3 s states at E^s ci - Ep = 11 .95 eV, 
and the one for the CI 4/ states at £'4y;ci - = 16.08 eV. We 
remark that the CI 3 p states give rise to small d and / contribu- 
tions to the DOS in K that are energetically far away from the 
E},d,K and £'4/^,k energy parameter, respectively. 

Fee Argon is the third investigated material. It also features 
a large band gap and enables the usage of a large MT radius 
that we set to Rmt = 3.15 «(>. We employ the experimental 
lattice constant of flexp = 9.93 ao, Gmax^Mx = 13 (-295 LAPW 
basis functions per atom), /^ax = 10, and 60 k points in the 
IBZ. Density of states and energy parameters for this material 
are shown in Figure |4] The energy parameters for the 3 s and 4/ 
states are beyond the energy range of the figure at E^s - Ep - 
-14.43 eV and £4/ - £f = 20.56 eV. Similarly to KCl, we 
discuss the band structure of Ar in detail in section l43] and show 
it in figure [141 
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Figure 4: Density of States and atomic energy parameters for Ar. 

To demonstrate the description of semicore states by the dif- 
ferent basis sets we also perform calculations on bcc vanadium 
with the experimental lattice constant Oexp - 5.73 qq. Here, we 
apply Rmt - 2.41 aq, Gmax = 4.5 Aq' (-145 LAPW basis func- 
tions per atom), /^ax = 10, and 190 k points in the IBZ. The 
energy parameters are at £4, - = -3.34 eV, E^p - Ep = 
-0.22 eV, E3d - = 0.05 eV, and £4/ - Ep ^ 2.05 eV. 

4. Linearization error 

In this section we quantify the linearization error for the dif- 
ferent basis sets. First, we define a basis representation error 
that quantifies the linearization error and examine its depen- 
dence on the energy parameters and MT sphere radii. Up to that 
point, we only consider non-selfconsistent calculations, where 
the ability of the basis to represent the exact single-particle 
wavefunctions in the spheres to a given effective potential is 
assessed. Then, we turn to self-consistent calculations. In par- 
ticular, we examine how total energies, lattice constants, and 
KS band gaps depend on the energy parameters and the MT 
sphere radii. We wiU see that the LO extension reduces these 
undesirable dependences to a great extent. 

4.1. MT basis representation 

In order to assess the quality of the basis representation in the 
MT spheres, we define the representation error 



A,{E) 



\\Ul - M/ll 



r^ [ui(r,E)-ui(r,E)Ydr 



(9) 



where M/(r, E) is a normalized solution of Eq. (|4|i for a given 
energy E and S/(r, E) is the best representation of M/(r, E) in 
terms of the radial functions of the given basis in the MT sphere, 
i.e, M/(r, E)-ui{r, E) is orthogonal to the basis. Thus, restrictions 
of the representation of m/ by &/ due to the matching to the plane 
waves in the IR at the MT boundary are not considered. The 
error A;(£') becomes zero if M/(r, E) can be represented point 
wise by the basis. This is the case if E is identical to an energy 




Figure 5: Representation en'or in the d channel of the MT sphere in fee Ce. 
There is a localized d core state at £ - E^j = -3.87 Htr, which is orthogonal to 
the space spanned by the basis sets. The inset shows a detailed view around the 
energy parameter in a logarithmic scale. 



parameter. On the other hand. A/ = 1 if ui(r, E) is orthogonal to 
the function space spanned by the available radial functions. 

In the following, we investigate the energy dependence of the 
representation error A/(£') for the conventional LAPW basis and 
its HDLO and HELO extensions. Based on the effective poten- 
tial from a self-consistent DFT calculation. Figure |5] shows the 
error for / = 2 in the vicinity of the energy parameter for the 5d 
states of bulk fee Ce. 

For an accurate representation of the electron density, which 
is the central quantity of DFT, the valence states must be de- 
scribed accurately. Therefore, the energy parameters are chosen 
in the valence band. The inset of Fig.|5]shows an energy interval 
of about ±0.3 Htr ^ ±8 eV around the Ce 5d energy parame- 
ter. In this energy range the linearization error of the conven- 
tional LAPW basis for the d channel of Ce amounts to maxi- 
mally 1.5 ■ 10"^. This error can be further reduced by adding 
LOs. While the HELOx 1 and HELOx2 basis sets reduce the 
maximal error to 3.3 ■ 10"^ and 1.1 ■ lO""*, respectively, the 
LAPW+HDLOxl exhibits a maximal error of only 3.8 ■ 10 a 
factor of 40 smaller than in the case of the conventional LAPW 
basis. We note that in the immediate vicinity of the energy pa- 
rameter, A;(£') scales as \E - /i/p for the conventional LAPW 
basis and as \E - Ei\^ for the HDLOxl basis. 

For energetically high-lying states the conventional LAPW 
basis quickly becomes inadequate, while LOs, especially the 
HELOs (cp. Fig. |5]l, can improve the basis substantially in a 
systematic and controllable manner. Methods that rely on the 
empty states such as the GW approximation [24] or RPA corre- 
lation functionals 125, 26] thus require an LAPW basis that is 
thoroughly converged with respect to LOs. 

For all basis sets we observe a sharp peak at -3.87 Htr 
{— -105 eV) where t^i{E) becomes 1, which signals a state that 
is orthogonal to the basis functions and therefore cannot be de- 
scribed by the basis. In practice, this is no problem as this is the 
Ad core state of Ce which is treated separately from the valence 
states. As it is completely confined to the MT sphere, it lies 
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Figure 7: Representation error at the fixed energy E^j - 0.1 Htr as a function 
Figure 6: Same as Fig. |5] for tlie p channel in bcc V. The semicore state at (he MT radius in fee Ce. The LAPW+HDLOxl basis suppresses the rep- 

-1.7 Htr is not completely confined in the MT sphere. Therefore, it is not or- resentation eiTor most efficiently. The error is below IQ-'. Also compare the 
thogonal to the space spanned by the basis sets and appears as shallow peaks logaiithmic plot in the inset 
with amplitudes well below 1 . 



outside the space spanned by the basis sets. On the other hand, 
the wave function of a semicore state extends significantly over 
the MT sphere. In particular, it does not vanish at the boundary 
and will therefore have a finite overlap with the valence basis. 
Figure |6] shows such a case in the p channel of bcc vanadium. 
The representation error features a shallow peak (or merely a 
shoulder) at -1.7 Htr (= -46 eV), which corresponds to the 
3p state. This precludes a treatment of the 3p state separate 
from the valence states, and we must employ LOs of / = 1 at 
the appropriate energy, especially if the LAPW+HDLOxl ba- 
sis is applied. It is evident that the rest of the diagram looks 
rather similar to the one in Fig. EJa), which demonstrates that 
the qualitative behavior of the representation error is essentially 
independent of the material and the angular momentum I. 

While the energy parameters can be adapted to the material 
in each iteration of the self-consistent-field run, the MT radius 
must be chosen before starting the calculation and then stays at 
this fixed value. One condition for the choice of the radius is 
that the MT spheres must not overlap. On the other hand, they 
should be chosen as large as possible in order to keep the re- 
ciprocal cutoff radius Gmax small. (The larger the spheres, the 
smaller the IR, and the less plane waves are necessary to rep- 
resent the wave functions there.) Finally, if one considers more 
than one chemical element, the sizes of the MT spheres relative 
to each other can be chosen according to tabulated atomic radii. 
This shows that there is no optimal choice of the radii. Ideally, 
the calculated results should be independent of the choice of the 
MT radii. 

In Fig.|2]we show the representation error foiSd states in Ce 
as a function of the MT radius for the different basis sets and at 
a fixed energy £ of 0. 1 Htr below the 5d energy parameter In 
the calculations we increased the cutoff Gmax while reducing the 
MT radius to keep the product Gmax^MX fixed to 13. The effec- 
tive potential was determined from self-consistent calculations 
for each MT radius. All other numerical parameters, as well as 
the method to choose the energy parameters, are identical to the 



previous Ce calculation. 

We observe that, independently of the basis chosen, the 
representation error becomes smaller as the MT radius is re- 
duced. This is related to the growing independence of M;(r, E) 
on the energy parameter E with decreasing radius observed in 
Fig. [Tib) and is most pronounced for the conventional LAPW 
basis, which has no additional LOs apart from those for the 
semicore states. From 3.15 aq to 2.4 aq the error decreases from 
1.3 ■ 10-^ to 2.5 ■ 10-4. In the case of the LAPW+HELO and 
LAPW-hHDLO basis sets, the dependence on the MT radius is 
strongly reduced. In particular, the LAPW+HDLOxl basis ex- 
hibits a very small representation error for this wide range of 
MT radii, nearly a factor 100 smaller than in the conventional 
LAPW basis. Such a stability with respect to non-convergence 
parameters like MT radii is clearly desired. 

4.2. Total energies and lattice constants 

Now we move to self-consistent calculations and address 
the questions whether the observations of the previous section 
translate to physical properties, such as the total energy Zstotai or 
the equilibrium lattice constant a after self-consistency in the 
electron density is achieved. Thus, we write the total energy 
£totai(a I {El], {^mt}) in a form where besides its dependence on 
the lattice constant a, the number of Bloch vectors, and the size 
of the basis set, which both have been chosen sufficiently large 
that the total energy can be considered converged with respect 
to these parameters, its dependence on the set of energy param- 
eters and MT radii becomes explicit. It is clear that the equilib- 
rium lattice constant a is obtained by minimization of the total 
energy with respect to variations of the lattice constant. The to- 
tal energy fitotai itself is dependent on the method-inherent pa- 
rameters. As in the previous section, we focus on the stability 
of the results with respect to variations of the energy parameters 
and the muffin-tin radii. 

In several respects, self-consistent calculations are more de- 
manding on the performance of basis sets than what we have 
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Figure 8: Total energy of fee Ce for the different basis sets from self-consistent- 
field calculations as a function of the energy parameter E1-2 relative to the 
Fernii energy. The total energy is plotted relative to the minimal energy ob- 
tained with the basis set including HDLOs. The inset shows a detailed view, 
reveahng that the LAPW+HDLOx 1 basis exhibits the smallest error by far. We 
also mark the energetic positions of the energy parameter as determined from 
two common approaches, which are desciibed in the text. In these calculations 
HDLOs and HELOs are only added in the d channel. 

considered so far. Firstly, we now account for the full non- 
spherical potential in the MT spheres, while we have restricted 
ourselves to the spherical potential before. Second, the IR is 
now taken into account as well, and the wave functions must be 
described accurately over the complete space. Third, the wave 
functions must also be described accurately over a sufficiently 
wide energy range to yield a precise valence electron density in 
each self-consistent-field iteration. 

In Fig. [8] we show the dependence of the calculated total en- 
ergy of fee Ce on the energy parameter £'/=2 for different basis 
sets, where, as an exception, the LAPW basis for the valence 
electrons is only extended by LOs of d character While keep- 
ing all other energy parameters fixed (to values taken from a 
previous self-consistent-field calculation), we vary Ei=2 in the 
range from -0.29 to 0.08 Htr -7.9 to 2.2 eV) relative to the 
Fermi level. As seen in the figure, the total energy calculated 
with the conventional LAPW basis set depends significantly on 
the choice of the energy parameter The marked energy param- 
eters determined from the two automatic approaches yield total 
energies that deviate substantially from each other and from the 
minimal total energy obtained with Ei=2-Ep - -0.16 Htr which 
lies about 1.3 eV below the lower valence band edge. 

When extending the LAPW basis with LOs, we observe that 
the total energy becomes lower (due to the increase of the vari- 
ational freedom enabled by the larger basis sets), but, more im- 
portantly, its dependence on £;=2 is strongly suppressed, most 
strongly again for the LAPW-i-HDLOxl set, which exhibits a 
variation that is two orders of magnitude smaller than for the 
conventional LAPW basis and, thus, negligible for all practi- 
cal purposes. Although the HELOx2 extension contains one 
set of LOs more than HDLOxl, the corresponding total energy 
exhibits a larger variation albeit still a very small one when 
compared to the conventional set. Of course, the calculations 




Figure 9: Ground-state total energy of fee Ce from self-consistent-field calcula- 
tions as a function of the MT radius for different basis sets. HDLOs and HELOs 
are added for / = 0, 1, 2, 3. The inset shows again that the LAPW+HDLOx 1 
basis exhibits the smallest error. 




Figure 10: Same as Fig.|9]for KCl. 



with the HELOs still depend on E^}^, i.e, the energy parame- 
ter of the HELOs. As akeady explained above, we employ for 
the HELOs energy parameters determined from Eq. dg). The 
HDLO accuracy could be realized with normal LOs if their en- 
ergy parameter Ef^ was chosen close to £'/=2- In fact, in the 
limit E^^ Ei^2 the two approaches are theoretically identi- 
cal. However, if the distance between the two energy param- 
eters gets too small, the near linear dependence of the basis 
might easily cause numerical problems, which is avoided by 
using HDLOs. 

Besides its dependence on the energy parameters, in Fig. |7] 
we already showed that the linearization error also strongly de- 
pends on the MT radii. The error decreases for smaller radii, 
but then the IR becomes larger, which entails a larger recipro- 
cal cutoff radius giving rise to an increase of the basis set size 
and thus higher computational demands. As we have seen, a 
more effective means to reduce or even eliminate the lineariza- 
tion error is the introduction of LOs, in particular the HDLOs. 
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Figure 1 1 : Equilibrium lattice constant as a function of the MT radius for fee 
Ce detemiined from Mumaghan fits to the total energies for 15 lattice constants 
in the range 0.938aexp to 0.966«[:xp- 
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Figure 12: Same as Fig. 1 11 I for KCl. Nine total-energy calculations for lattice 
constants in the range 0.958flexp to 0.974aexp were employed for the fit. 



We observe the same behavior for the total energy and the 
equilibrium lattice constant, for which the total energy assumes 
a minimum. The latter is calculated from a series of total- 
energy calculations for different lattice constants together with 
a Murnaghan fit 12711 . We show results for Ce and KCl. To 
avoid overcompleteness of the basis at small MT radii, we nei- 
ther include HDLOs nor HELOs in the s and p channels of K, 
where we already have LOs to describe the semicore states. 

In Figs. |9] and [TOl the total energy is shown as a function of 
the MT radius. In fact, the diagrams look qualitatively very 
similar to the MT representation error shown in Figs. 2]and|7] 
respectively. The variation in the total energy is between one 
and two orders of magnitude smaller with the HDLO extension 
than without for both compounds. 

Even without the LOs, the total energy is accurate up to few 
mHtr. However, the equilibrium lattice constant determined 
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Figure 13: Band structure of KCl calculated with the conventional LAPW and 
the LAPW+HDLOx 1 basis set. On the given energy scale the band structures 
calculated with the HELO extensions are indistinguishable from the one with 
HDLOs. 



from the minimum of the total energy depends on small en- 
ergy differences. It is thus particularly susceptible to inaccu- 
racies in the total energies. As shown in Fig. fTTIandfT^] the 
variation of the lattice constant a is on the order of 1 percent 
for Ce (flexp = 9.05 ao, alapw = 8.55 to 8.69 ao) and KCl 
(flexp = 11.89 flo, flLAPW = 1 1 .49 to 1 1 .54 flo) and thus in the 
same ballpark as the supposed accuracy of the LDA functional. 
After addition of one set of HDLOs the variation of a is strongly 
suppressed to less than 0.1% (Ce: alapw+holoxi - 8.539 to 
8.544 flo, KCl: flLAPW+HDLOxi = 11.487 to 11.488 flo). 

We note that each additional set of LOs (for / = 0, . . ., 3) in- 
creases the basis-set size by only 16 functions per atom. Thus, 
adding LOs is computationally much more efficient than reduc- 
ing the MT radius, which would make many more augmented 
plane waves necessary to enable an adequate description of the 
wave functions in the interstitial region. 

4.3. KS band gap 

The band gap is a fundamental quantity of semiconductors 
and insulators and is given by the energy difference of the low- 
est unoccupied and the highest occupied state. Its calculated 
value also depends on the ability of the basis set to describe 
both, the valence and the conduction bands with sufficient ac- 
curacy. Since the energy parameters are located typically in the 
energy range of the valence bands (see Sec. 14.2b . one might ex- 
pect deviations of the conduction states from the true KS eigen- 
states having an adverse effect on the accuracy of the band-gap 
value. We demonstrate for KCl and Ar that the inaccuracies of 
the KS gap due to the linearization error can indeed be sizable 
unless LOs are used to ensure a proper description in the MT 
spheres. With the LAPW basis set, the KS band gap is calcu- 
lated to be larger than with the extended basis sets, leading to 
an underestimation of the KS band gap error when compared to 
experimental values. 
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Figure 14: Same as Fig. 1 1 3 I for Ar. 



In Figs. [13] and [H the KS band structures of rock-salt KCl 
and crystalline Ar, respectively, are shown. Both exhibit a 
direct band gap at T. We show band structures calculated 
with the conventional LAPW and the LAPW+HDLOxl ba- 
sis (results obtained by the HELO extensions lie on top of the 
LAPW+HDLOxl ones on the scale of the diagrams). As ex- 
pected, the occupied states, being close to the energy parame- 
ters, are already well described in the conventional LAPW ba- 
sis and the occupied bands are indistinguishable on the energy 
scale of the diagrams. However, there are clear deviations in the 
unoccupied states. In particular, the lowest unoccupied band, 
which is of 4s character, shifts down by 0.19 eV for KCl and 
even by 1 .87 eV for Ar upon introducing the LOs; the band dis- 
persion is also affected noticeably. The downward shift of the 
lowest unoccupied band entails, of course, a reduction of the 
KS band gaps, which amounts to 4% for KCl and 19% for Ar. 
As a result, the linearization error has an even more pronounced 
effect on the band gap than it has on the total energy and equi- 
librium lattice constant. Higher-lying unoccupied states will be 
affected even more strongly so that a thorough convergence of 
the basis set is mandatory in the calculation of quantities that 
involve the unoccupied states in a perturbative way such as re- 
sponse quantities used, e.g., in the optimized-effective-potential 
method and the GW approximation Ii20n23il . 

5. Basis-set size 

In this section, we compare for the example of fee Ce the 
speed of convergence of the total energy with respect to the 
cutoff radius of reciprocal lattice vectors achieved with the dif- 
ferent basis sets. In Fig.[T5]the convergence of the total energy 
is shown as a function of Gmax^MX with Rmt = 3.14 aq- In ad- 
dition to the conventional LAPW basis and those extended by 
HDLOs and HELOs, we also consider the convergence for the 
APW+lo basis and its HELOxl extension. In contrast to the 
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LAPW basis, the APW+lo approach 112 111 augments the inter- 



Figure 15: Convergence of the total energy for fee Ce with respect to Gmax^MT 
withfijyiT = 3.14«o- The total energy is given relative to the minimal energy ob- 
tained with the LAPW+HDLOx 1 basis. The upper abscissa of the figure gives 
the size of the conventional LAPW basis averaged over all k points. For basis 
sets with one set of LOs (dashed fines, LAPW+HELOxl, LAPW+HDLOxl, 
APW+lo) the total number of basis functions is increased by 16, for those with 
two sets (dotted lines. LAPW+HELOx2, APW+lo+HELOxl) it is increased 
by 32. The best performance in terms of rate of convergence and accuracy of 
the converged value is achieved with the LAPW+HDLOxl basis. 



stitial plane waves solely by the functions M/(r, £'/)F/,„(f ), while 
the energy derivatives M/(r, £/) up to a given angular momentum 
(here, / = 3) are treated as LOs (for which we have adopted 
the common lower case notation "lo"). For higher I we de- 
part from the ori gina l APW+lo approach and follow the idea 
of Madsen et al. 112 811 to use both radial functions for augmen- 
tation as in LAPW. As a consequence of the altered lineariza- 
tion, the APW+lo basis is only continuous in value but shows 
a discontinuity in the slope at the MT sphere boundary, i.e., the 
basis functions exhibit a kink there. The less stringent match- 
ing at the MT sphere boundary leads to a more flexible basis set, 
which manifests in a faster convergence of the total energy with 
respect to Gmax^Mi in comparison to the conventional LAPW 
basis. This effect is clearly visible in Fig. [15] While the conven- 
tional LAPW basis does not yield a converged total energy even 
with Gniax^MT = 13, the APW+lo curve reaches its converged 
value already for Gmax^MX ~ 10. We note that the LAPW basis 
finally converges to the same total energy as the APW+lo ap- 
proach, but at a considerably larger basis. This is in agreement 
with earlier investigations on fee Ce performed by Sjostedt et 
al. ||2lll . 

However, the converged values of both LAPW and APW+lo 
calculations deviate substantially from the lower and more 
accurate total energies obtained with the basis sets that are 
extended with additional LOs. Furthermore, the APW+lo 
curve should rather be compared with LAPW+HELOxl or 
LAPW+HDLOxl, which have the same basis-set size for a 
given Gmax- The latter basis reaches the variationally lowest 
total energy among all basis sets. We note that adding further 
HDLOs (i.e., the third and higher derivatives) hardly changes 
the converged value. We can also observe that the correspond- 
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ing calculations converge more rapidly than those employing 
the conventional LAPW basis. The reason for that is the in- 
creased flexibility brought about by the additional LOs. This 
increased flexibility plays a similar role as allowing for a kink 
in the APW+lo basis functions. So, it acts toward a decou- 
pling of the two regions of space, allowing to some extent a 
separate variational search for the eigenstates in the two re- 
gions. The same general observations can be made when com- 
paring the APW+lo+HELOxl and the LAPW+HELOxl ba- 
sis sets. The two curves converge to the same total-energy 
value, which is, however, still higher than that obtained with 
the LAPW+HDLOxl basis. It should be pointed out that ev- 
ery kinkless wave function that is given as a linear combination 
of plane waves in the IR and the M/(r, £/) and M/(r, £;) in the 
MT spheres is representable by LAPW basis functions since 
the wave function has to feature the MT boundary matching 
that is naturally obtained with the LAPW basis. Hence, when- 
ever an APW+lo calculation yields a smaller total energy than 
the corresponding LAPW calculation, then the wave functions 
necessarily exhibit a kink at the MT sphere boundaries, which 
from a physical point of view might be unsatisfactory. Loosely 
speaking, the variational search for the wave functions reached 
a lower total energy by accepting a kink with a singular kinetic 
energy density at the sphere boundaries. 

6. Conclusion 

In foresight of new challenges treating solids with greater 
electronic and chemical complexity with more sophisticated 
functionals, we have explored in this work the capability of the 
LAPW basis set to deal with these challenges and evaluated 
extensions of the basis set by two types of local orbitals that 
basically offer the potential to provide the density functional 
answer to the problem at hand with very high precision. 

In detail we have analyzed the effects of the linearization er- 
ror that arises in the LAPW basis due to the restriction to only 
two radial functions, «;(r, £/) and «/(r, £/), per / quantum num- 
ber in the MT spheres. While the LAPW basis is established 
to be very accurate for a wide range of materials and material 
properties, there are cases where the linearization error c auses 
the results to depend appreciably on numerical parameters that 
are inherent to the FLAPW method and are not convergence 
parameters, i.e., the energy parameters and the MT radii. Al- 
though choosing small MT radii reduces the linearization error, 
this extends the IR and thus entails a larger reciprocal cutoff" ra- 
dius leading to many additional basis functions. A much more 
effective way is to add LOs, for example HELOs, which are de- 
fined with higher energy parameters. HELOs are often used to 
improve the description of unoccupied states, but can also be 
employed to greatly reduce the linearization error for the oc- 
cupied states. An even more efficient way to eliminate the lin- 
earization error for the occupied states, and thereby make the 
results insensitive to variations of the energy parameters and 
MT radii, is to extend the basis with HDLOs, which are con- 
structed from the second energy derivatives M/(r). 

The addition of LOs allows to employ large MT spheres, 
which keeps the reciprocal cutoff" radius small and thus keeps 



the basis-set size to a minimum. Using large spheres also en- 
sures that the orthogonality of the basis functions for the va- 
lence states to the core states is maximized, which avoids ghost 
bands to appear in the valence and conduction band structure. 

Similar to the APW-i-lo approach, extending the conventional 
LAPW basis set by LOs leads to a more ffexible basis set result- 
ing in a faster convergence with basis-set size. In cases where 
the linearization error is large, the LAPW+HDLOxl basis ex- 
hibits higher precision at equal basis set sizes in comparison to 
APW+lo. 
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